#---------------------------------------#
#---------------------------------------#
# Data Analysis Script
# Article: Partisan responses to democracy promotion – Estimating the causal effect of a civic information portal
# Journal: World Development
# Date: 2020-02-10
#---------------------------------------#
#---------------------------------------#

#---------------------------------------#
# Packages
#---------------------------------------#

library(zoo)
library(mfx)
library(ri)
require(memisc)
library(pander)
library(stargazer)
require(gridExtra)
library(jtools)
library(sjPlot)
library(officer)
library(flextable)

# Options
options(scipen=999)


#---------------------------------------#
#Data
#---------------------------------------#

path <- c("[path]")
setwd(path)
load("Mzalendo_RCT_data.RData")

# Path for graophs
graphs.path <- c("[path]")
setwd(graphs.path)


#---------------------------------------#
# Coding
#---------------------------------------#

# Subset
# subset <- c("Non-Part-in-Govt","Non-Part-in-NGovt","Undecl-in-Govt","Undecl-in-Non-Govt")  # DELETE
subset <- c("Oppos-in-Govt","Govt-in-Govt","Govt-in-Non-Govt","Oppos-in-Non-Govt")


#---------------------------------------#
# Balance & Attrition
#---------------------------------------#

#---------------------------------------#
# Input for Figure 2. Flow Diagram of Subjects Over the Course of the Experiment

# Wave1
raw$wave1.start <- ifelse(is.na(raw$s1_q1)==T,F,T)  # Includes incompletes, that started, but did not finnish
table(raw$wave1.start)
raw$wave1.completed <- ifelse(raw$email=="",F,T)
table(raw$wave1.completed)

# Treatment assignment
table(raw$treatment[raw$wave1.completed])

# Wave2
raw$wave2.start <- ifelse(is.na(raw$s2_q4)==T,F,T) #includes incompletes, that started, but didn't finnish
table(raw$wave2.start)
tapply(raw$wave2.start, raw$treatment, sum, na.rm=T)
raw$wave2.completed <- ifelse(is.na(raw$s2_q15)==T,F,T)
raw$wave2.completed[raw$s2_q15==""] <- F
table(raw$wave2.completed)
tapply(raw$wave2.completed, raw$treatment, sum, na.rm=T)
sum(tapply(raw$wave2.completed, raw$treatment, sum, na.rm=T))


#---------------------------------------#
# Appendix 2: Covariate Balance
# Figure 5. Covariate Balance - F-Test Simulations To Assess The Null Hypothesis That The Covariates Predict Random Assignment (Z) No Better Than Would Be Expected By Chance.

set.seed(1234567)
vars <- c("treatment","time_started.s1","age","female","uni","ses1","pro.govt")
data.sim <- data[complete.cases(data[vars])==T,vars]
for (i in vars) {
	data.sim[,i] <- as.integer(data.sim[,i])
}  
numiter <- 10000
Z <-  data.sim$treatment
covs <- as.matrix(data.sim[,vars[-1]])
perms <- genperms(Z,maxiter=numiter)
# F-statistic from actual data.sim
Fstat <- summary(lm(Z~covs))$fstatistic[1]
# Simulate
Fstatstore <- rep(NA,numiter)
for (i in 1:numiter) {
	Fstatstore[i] <- summary(lm(perms[,i]~covs))$fstatistic[1]   # F-statistic under the null of random assignment of Z
}
# Plot
hist(Fstatstore,freq=TRUE,xlab="Estimated F-stat",
		 main=paste("Distribution of the F-statistic\np = ",
		 					 round(sum(Fstatstore >= Fstat)/numiter,3),", SE = ",
		 					 round(sd(Fstatstore),3),sep=""), breaks=numiter,lwd=1)
lines(x=c(Fstat,Fstat),y=c(-1000,numiter*2),lwd=2,col="red",lty=2)


#---------------------------------------#
# Appendix 8: Additional Statistical Material

# Table 5. Covariate Balance per Main Partisanship Category – Predictors of Treatment Assignment, Linear Models (OLS).
summary(lm(treatment ~ age + female + uni + ses1 + pro.govt.full, data))
balance.w2.govgov <- lm(treatment ~ age + female + uni + ses1, 
												data[data$pro.govt.full %in% "Govt-in-Govt",])
balance.w2.govopp <- lm(treatment ~ age + female + uni + ses1, 
												data[data$pro.govt.full %in% "Govt-in-Non-Govt",])
balance.w2.oppgov <- lm(treatment ~ age + female + uni + ses1, 
												data[data$pro.govt.full %in% "Oppos-in-Govt",])
balance.w2.oppopp <- lm(treatment ~ age + female + uni + ses1, 
												data[data$pro.govt.full %in% "Oppos-in-Non-Govt",])
coef.names <- c("(Intercept)","Age","Female","University","Socio-Economic")
names(balance.w2.govgov$coefficients) <- coef.names
names(balance.w2.govopp$coefficients) <- coef.names
names(balance.w2.oppgov$coefficients) <- coef.names
names(balance.w2.oppopp$coefficients) <- coef.names
balance.w2.subsets <- mtable('Gov-Gov' = balance.w2.govgov,
														 'Gov-Opp' = balance.w2.govopp,
														 'Opp-Gov' = balance.w2.oppgov,
														 'Opp-Opp' = balance.w2.oppopp,
														 summary.stats = c('N'))
show_html(balance.w2.subsets)
write.mtable(balance.w2.subsets,file="balance_subsets.csv", colsep=",")


# Table 6. Prediction of Completion of the Endline Survey, Linear Models (OLS).
table(raw$wave2.completed[raw$wave1.completed],
				 raw$treatment[raw$wave1.completed])
1-tapply(raw$wave2.completed[raw$wave1.completed],
				 raw$treatment[raw$wave1.completed], mean)
fit.wave2.start <- lm(wave2.start ~ treatment, data=raw[raw$wave1.completed,])
fit.wave2.completed <- lm(wave2.completed ~ treatment, data=raw[raw$wave1.completed,])

coef.names <- c("(Intercept)","Treatment")
names(fit.wave2.start$coefficients) <- coef.names
names(fit.wave2.completed$coefficients) <- coef.names
attrition <- mtable('fit.wave2.start' = fit.wave2.start,
														 'fit.wave2.completed' = fit.wave2.completed,
														 summary.stats = c('N'))
write.mtable(attrition,file="attrition.csv", colsep=",")


#---------------------------------------#
# Descriptive Statistics
#---------------------------------------#

#---------------------------------------#
# Table 1. Descriptive Statistics for the Subject Pool.
stargazer(data[,c("age","female","ses1","uni","ext.efficacy.w2","s2_q7")], 
					type = "html", title="Descriptive statistics", digits=2,
					covariate.labels=c("Age","Female","Socio-Economic (wealth index)","University (attended some)",
														 "Extrnal Efficacy Score (1-5)","Likelihood of Voting (1-4)"),
					header=F, out="descriptive.txt")


#---------------------------------------#
# Figure 4. Partisanship Categories – All Observations. 
set_theme(
	base = theme_blank(),
	axis.title.size = 1,
	axis.textsize.y = 1,
	axis.textsize.x = 1,
	legend.size = .6,
	geom.label.size = 4,
	axis.angle.x = 90
)
plot_frq(data$pro.govt.full,
				 axis.title = "")


#---------------------------------------#
#  Figure. External Efficacy - Descriptive
likert.function <- function (vars) {
	factor(vars, levels=1:5, labels=c("Strongly Disagree", "Disagree",
																		"Neither", "Agree", "Strongly Agree"))
}
library(likert)
#subset and rename
items <- data[which(data$s2_q15!=""),c("ext.eff1.w1","ext.eff2.w1", "ext.eff3.w1")] #Qs of interest
for(i in 1:ncol(items)) {
	items[,i] = likert.function(items[,i])
}
names(items) = c("1: Care what people think", 
								 "2: Interest only in votes (inverse)", 
								 "3: Responsive to needs")
plot(likert(items), group.order=names(items), low.color = "#D8B365", high.color = "#5AB4AC",
		 neutral.color = "grey90", neutral.color.ramp = "white") +
	theme(legend.text=element_text(size=13), legend.title=element_blank())



#---------------------------------------#
# Figure. Participation index likert
likert.function <- function (vars) {
	factor(vars, levels=1:4, labels=c("Very Unlikely", "Somewhat Unlikely",
																		"Somewhat Likely", "Very Likely"))
}
library(likert)
#subset and rename
items <- data[which(data$s2_q15!=""),c("s1_q6","s1_q7")] #Qs of interest
for(i in 1:ncol(items)) {
	items[,i] = likert.function(items[,i])
}
names(items) = c("Contact MP (Q6)", 
								 "Vote in Elections (Q7)")
plot(likert(items), low.color = "#D8B365", high.color = "#5AB4AC",
		 neutral.color = "grey90", neutral.color.ramp = "white") +
	theme(legend.text=element_text(size=13), legend.title=element_blank())


#---------------------------------------#
# Results section: External efficacy
#---------------------------------------#

# PAP model
levels(data$pro.govt.raw)
fit.ext.int1 <- lm(ext.efficacy.w2 ~ treatment + treatment:pro.govt.raw + age + female + ses1 + uni + pro.govt.raw + ext.efficacy.w1, data)
# Const. dummy
fit.ext.int2 <- lm(ext.efficacy.w2 ~ treatment + treatment:pro.govt.in.govt.na + age + female + ses1 + uni + pro.govt.in.govt.na + ext.efficacy.w1, data)
# Const. full
levels(data$pro.govt.full)
fit.ext.int3 <- lm(ext.efficacy.w2 ~ treatment + treatment:pro.govt.full + age + female + ses1 + uni + pro.govt.full + ext.efficacy.w1, data)

# Rename coef
coef_rename <- function(model_name) {
	coefs <- names(model_name$coefficients)
	coefs <- gsub('pro.govt.raw', '', coefs)
	coefs <- gsub('pro.govt.in.govt.naTRUE', 'Govt-in-Govt', coefs)
	coefs <- gsub('pro.govt.full', '', coefs)
	coefs <- gsub('treatmentTRUE', 'Treatment', coefs)
	coefs <- gsub(':', '*', coefs)
	coefs <- ifelse(coefs=='Treatment', 'Treatment (ref.cat)', coefs)
	coefs <- ifelse(coefs=='age', 'Age', coefs)
	coefs <- ifelse(coefs=='femaleTRUE', 'Female', coefs)
	coefs <- ifelse(coefs=='ses1', 'Socio-Econ', coefs)
	coefs <- ifelse(coefs=='uniTRUE', 'University', coefs)
	coefs <- ifelse(coefs=='ext.efficacy.w1', 'Baseline efficacy', coefs)
	return(coefs)
}
names(fit.ext.int1$coefficients) <- coef_rename(fit.ext.int1)
names(fit.ext.int2$coefficients) <- coef_rename(fit.ext.int2)
names(fit.ext.int3$coefficients) <- coef_rename(fit.ext.int3)

export_summs(fit.ext.int1, fit.ext.int2, fit.ext.int3,
						 omit.coefs = '(Intercept)',
						 scale = F, to.file = "docx", file.name = "ext_eff_models.docx")
export_summs(fit.ext.int3,
						 scale = F, to.file = "pdf", file.name = "ext_eff_models.pdf")

# Visual
plot_summs(fit.ext.int3,
					 coefs = c(names(fit.ext.int3$coefficients)[grep('Treat*',names(fit.ext.int3$coefficients))]),
					 scale = F)

# Substantive interpretation
(-0.046813+0.222339)/2.48


#---------------------------------------#
# Explore partial winners in more detail
# [In text references to Statistical Appendix]
# Table 9. Treatment Effect on Diffuse and Constituency-Specific External Efficacy, Linear Model (OLS).
data$ext.efficacy.diffuse.w1 <- (rowSums(data[,c("ext.eff1.w1","ext.eff2.w1")])/2)
data$ext.efficacy.diffuse.w2 <- (rowSums(data[,c("ext.eff1.w2","ext.eff2.w2")])/2)

# Const. full
diffuse.int <- lm(ext.efficacy.diffuse.w2 ~ treatment + treatment:pro.govt.full + age + female + ses1 + uni + pro.govt.full + ext.efficacy.diffuse.w1, data)
specific.int <- lm(ext.eff3.w2 ~ treatment + treatment:pro.govt.full + age + female + ses1 + uni + pro.govt.full + ext.eff3.w1, data)

# Rename coef
coef_rename <- function(model_name) {
	coefs <- names(model_name$coefficients)
	coefs <- gsub('pro.govt.raw', '', coefs)
	coefs <- gsub('pro.govt.in.govt.naTRUE', 'Govt-in-Govt', coefs)
	coefs <- gsub('pro.govt.full', '', coefs)
	coefs <- gsub('treatmentTRUE', 'Treatment', coefs)
	coefs <- gsub(':', '*', coefs)
	coefs <- ifelse(coefs=='Treatment', 'Treatment (ref.cat)', coefs)
	coefs <- ifelse(coefs=='age', 'Age', coefs)
	coefs <- ifelse(coefs=='femaleTRUE', 'Female', coefs)
	coefs <- ifelse(coefs=='ses1', 'Socio-Econ', coefs)
	coefs <- ifelse(coefs=='uniTRUE', 'University', coefs)
	coefs <- ifelse(coefs=='ext.efficacy.w1', 'Baseline efficacy', coefs)
	coefs <- ifelse(coefs=='ext.efficacy.diffuse.w1', 'Baseline efficacy', coefs)
	coefs <- ifelse(coefs=='ext.eff3.w1', 'Baseline efficacy', coefs)
	return(coefs)
}
names(diffuse.int$coefficients) <- coef_rename(diffuse.int)
names(specific.int$coefficients) <- coef_rename(specific.int)

# New packages
export_summs(diffuse.int, specific.int,
						 coefs = c(names(diffuse.int$coefficients)[grep('Treat*',names(diffuse.int$coefficients))]),
						 scale = F, to.file = "pdf", file.name = "ext_eff_models.pdf")
export_summs(diffuse.int, specific.int,
						 coefs = c(names(diffuse.int$coefficients)[grep('Treat*',names(diffuse.int$coefficients))]),
						 scale = F, to.file = "xlsx", file.name = "diffuse_ext_eff_models.xlsx")


#---------------------------------------#
# Results section: Turnout Intention
#---------------------------------------#

#---------------------------------------#
# Table 3. Treatment Effect on Turnout Intention per Partisan Category, Separate Linear Models (OLS).

# PAP model
fit.partic.int1 <- lm(particip.w2 ~ treatment + treatment:pro.govt.raw + age + female + ses1 + uni + pro.govt.raw + particip.w1, data)
# Const. dummy
fit.partic.int2 <- lm(particip.w2 ~ treatment + treatment:pro.govt.in.govt.na + age + female + ses1 + uni + pro.govt.in.govt.na + particip.w1, data)
# Const. full
fit.partic.int3 <- lm(particip.w2 ~ treatment + treatment:pro.govt.full + age + female + ses1 + uni + pro.govt.full + particip.w1, data)

# Rename coef
coef_rename <- function(model_name) {
	coefs <- names(model_name$coefficients)
	coefs <- gsub('pro.govt.raw', '', coefs)
	coefs <- gsub('pro.govt.in.govt.naTRUE', 'Govt-in-Govt', coefs)
	coefs <- gsub('pro.govt.full', '', coefs)
	coefs <- gsub('treatmentTRUE', 'Treatment', coefs)
	coefs <- gsub(':', '*', coefs)
	coefs <- ifelse(coefs=='Treatment', 'Treatment (ref.cat)', coefs)
	coefs <- ifelse(coefs=='age', 'Age', coefs)
	coefs <- ifelse(coefs=='femaleTRUE', 'Female', coefs)
	coefs <- ifelse(coefs=='ses1', 'Socio-Econ', coefs)
	coefs <- ifelse(coefs=='uniTRUE', 'University', coefs)
	coefs <- ifelse(coefs=='ext.efficacy.w1', 'Baseline efficacy', coefs)
	coefs <- ifelse(coefs=='particip.w1', 'Baseline participation', coefs)
	return(coefs)
}
names(fit.partic.int1$coefficients) <- coef_rename(fit.partic.int1)
names(fit.partic.int2$coefficients) <- coef_rename(fit.partic.int2)
names(fit.partic.int3$coefficients) <- coef_rename(fit.partic.int3)

# New packages
export_summs(fit.partic.int1, fit.partic.int2, fit.partic.int3,
						 omit.coefs = '(Intercept)',
						 scale = F, to.file = "pdf", file.name = "partic_models.pdf")

# Visual
plot_summs(fit.partic.int3,
					 coefs = c(names(fit.partic.int3$coefficients)[grep('Treat*',names(fit.partic.int3$coefficients))]),
					 scale = F)

# COMBO Visual
plot_summs(fit.ext.int3, fit.partic.int3,
					 coefs = c(names(fit.ext.int3$coefficients)[grep('Treat*',names(fit.ext.int3$coefficients))]),
					 scale = F, model.names = c("Efficacy", "Participation"))

# Substantive interpretation
(-0.046813+0.222339)/3.57
# tapply(data$s1_q6, data$pro.govt.full, mean)

# The FDR-adjustment is done for all six hypotheses,
# three (winner, loser, partial winner) for both outcomes of interest (efficacy and participation)
# Create joint table
fdr.tab <- data.frame(rbind(
	summary(fit.ext.int3)$coefficients[c("Treatment (ref.cat)", "Treatment*Govt-in-Govt"),1:4],
	summary(fit.partic.int3)$coefficients[c("Treatment (ref.cat)", "Treatment*Govt-in-Govt"),1:4]))
rownames(fdr.tab) <- c("Efficacy: Opp-in-Gov","Efficacy: Gov-in-Gov",
											 "Participation: Opp-in-Gov","Participation: Gov-in-Gov")
fdr.tab$adjusted.p <- p.adjust(fdr.tab[,4], method="fdr")
colnames(ext.eff.tab) <- c("N","Estimate","SE","p-value","FDR-adj p")
ext.eff.tab <- round(ext.eff.tab, 4)
row.names(ext.eff.tab) <- c("Winner (Gov-Gov)","Loser (Opp-Gov)","Partial Winner (Opp-Opp & Gov-Opp)")


#---------------------------------------#
# Figure 3. Turnout Intention Likert Scale by Treatment and Control for the Two Partial Loser Partisans.
# Visualize Likert per subset
set_theme(theme = "blank",
					axis.title.size = .9,
					axis.textsize = .7,
					legend.size = 1,
					legend.title.size = 1,
					geom.label.size = 2,
					title.size = 1, 
					legend.inside = T,
					legend.pos = c(.3, .7),  # relative legend position inside plot
					legend.just = c(.1, .8))

subset.df <- data[data$pro.govt.full %in% "Govt-in-Non-Govt",]
s2_q7.graph.govopp <- sjp.xtab(subset.df$s2_q7, subset.df$treatment, show.total=F, axis.title="",
															 legend.labels=c("Control","Treatment"), legend.title="",
															 axis.labels=c("Very unlikely","Somewhat unlikely",
															 							"Somewhat likely","Very likely"),
															 title="Gov-Opp")
subset.df <- data[data$pro.govt.full %in% "Oppos-in-Non-Govt",]
s2_q7.graph.oppopp <- sjp.xtab(subset.df$s2_q7, subset.df$treatment, show.total=F, axis.title="",
															 show.legend=F, axis.labels=c("Very unlikely","Somewhat unlikely",
															 														 "Somewhat likely","Very likely"),
															 title="Opp-Opp")

grid.arrange(s2_q7.graph.govopp, s2_q7.graph.oppopp, nrow=1, ncol=2)


#---------------------------------------#
# Table 4. Intent to Contact MP in own Constituency, OLS regressions for Separate Partisan Categories.

# PAP model
fit.contactmp.int1 <- lm(s2_q6 ~ treatment + treatment:pro.govt.raw + age + female + ses1 + uni + pro.govt.raw + s1_q6, data)
# Const. dummy
fit.contactmp.int2 <- lm(s2_q6 ~ treatment + treatment:pro.govt.in.govt.na + age + female + ses1 + uni + pro.govt.in.govt.na + s1_q6, data)
# Const. full
fit.contactmp.int3 <- lm(s2_q6 ~ treatment + treatment:pro.govt.full + age + female + ses1 + uni + pro.govt.full + s1_q6, data)

mtable(
	'fit.contactmp.int3' = fit.contactmp.int3,
	summary.stats = c('N'))

# Rename coef
coef_rename <- function(model_name) {
	coefs <- names(model_name$coefficients)
	coefs <- gsub('pro.govt.raw', '', coefs)
	coefs <- gsub('pro.govt.in.govt.naTRUE', 'Govt-in-Govt', coefs)
	coefs <- gsub('pro.govt.full', '', coefs)
	coefs <- gsub('treatmentTRUE', 'Treatment', coefs)
	coefs <- gsub(':', '*', coefs)
	coefs <- ifelse(coefs=='Treatment', 'Treatment (ref.cat)', coefs)
	coefs <- ifelse(coefs=='age', 'Age', coefs)
	coefs <- ifelse(coefs=='femaleTRUE', 'Female', coefs)
	coefs <- ifelse(coefs=='ses1', 'Socio-Econ', coefs)
	coefs <- ifelse(coefs=='uniTRUE', 'University', coefs)
	coefs <- ifelse(coefs=='ext.efficacy.w1', 'Baseline efficacy', coefs)
	coefs <- ifelse(coefs=='particip.w1', 'Baseline participation', coefs)
	coefs <- ifelse(coefs=='s1_q6', 'Baseline participation', coefs)
	return(coefs)
}
names(fit.contactmp.int3$coefficients) <- coef_rename(fit.contactmp.int3)

# New packages
export_summs(fit.contactmp.int3, 
						 omit.coefs = '(Intercept)',
						 scale = F, to.file = "pdf", file.name = "contactmp_models.pdf")

# Visual
plot_summs(fit.contactmp.int3,
					 coefs = c(names(fit.contactmp.int3$coefficients)[grep('Treat*',names(fit.contactmp.int3$coefficients))]),
					 scale = F)

# Substantive interpretation
(-0.046813+0.222339)/3.57


#---------------------------------------#
# Table 10. Turnout Intention for Partial Loser, Linear Model (OLS).

turnout.govopp <- lm(s2_q7 ~ treatment + s1_q7 + age + 
					 						female + ses1 + uni, data[data$pro.govt.full %in% "Govt-in-Non-Govt",])
turnout.oppopp <- lm(s2_q7 ~ treatment + s1_q7 + age + 
										 	female + ses1 + uni, data[data$pro.govt.full %in% "Oppos-in-Non-Govt",])
# For table
coef.names <- c("(Intercept)","Treatment", "Turnout Intent (baseline)","Age","Female",
								"Socio-Economic","University")
names(turnout.govopp$coefficients) <- coef.names
names(turnout.oppopp$coefficients) <- coef.names
turnout.partlosers <- mtable('turnout.govopp' = turnout.govopp,
			 'turnout.oppopp' = turnout.oppopp,
			 summary.stats = c('N'))


#---------------------------------------#
# Table 11. Treatment Effect on External Efficacy Index per Other Partisan Categories, Sepa-rate Linear Models (OLS).

#External Efficacy
ext.eff.npgov <- lm(ext.efficacy.w2 ~ treatment + ext.efficacy.w1 + age + 
										 	female + ses1 + uni, data[data$pro.govt.full %in% "Non-Part-in-Govt",])
ext.eff.npopp <- lm(ext.efficacy.w2 ~ treatment + ext.efficacy.w1 + age + 
											female + ses1 + uni, data[data$pro.govt.full %in% "Non-Part-in-NGovt",])
ext.eff.udgov <- lm(ext.efficacy.w2 ~ treatment + ext.efficacy.w1 + age + 
											female + ses1 + uni, data[data$pro.govt.full %in% "Undecl-in-Govt",])
ext.eff.udopp <- lm(ext.efficacy.w2 ~ treatment + ext.efficacy.w1 + age + 
											female + ses1 + uni, data[data$pro.govt.full %in% "Undecl-in-Non-Govt",])
coef.names <- c("(Intercept)","Treatment", "Ext.Efficacy (baseline)","Age","Female",
								"Socio-Economic","University")
names(ext.eff.npgov$coefficients) <- coef.names
names(ext.eff.npopp$coefficients) <- coef.names
names(ext.eff.udgov$coefficients) <- coef.names
names(ext.eff.udopp$coefficients) <- coef.names
# Create joint table
ext.eff.other.part <- mtable(ext.eff.npgov, ext.eff.npopp, ext.eff.udgov, ext.eff.udopp,
			 summary.stats=c("N"))


#---------------------------------------#
# Table 12. Treatment Effect on Turnout Intention per Other Partisan Categories, Separate Linear Models (OLS)

turnout.npgov <- lm(s2_q7 ~ treatment + s1_q7 + age + 
											female + ses1 + uni, data[data$pro.govt.full %in% "Non-Part-in-Govt",])
turnout.npopp <- lm(s2_q7 ~ treatment + s1_q7 + age + 
											female + ses1 + uni, data[data$pro.govt.full %in% "Non-Part-in-NGovt",])
turnout.udgov <- lm(s2_q7 ~ treatment + s1_q7 + age + 
											female + ses1 + uni, data[data$pro.govt.full %in% "Undecl-in-Govt",])
turnout.udopp <- lm(s2_q7 ~ treatment + s1_q7 + age + 
											female + ses1 + uni, data[data$pro.govt.full %in% "Undecl-in-Non-Govt",])
coef.names <- c("(Intercept)","Treatment", "turnouticacy (baseline)","Age","Female",
								"Socio-Economic","University")
names(turnout.npgov$coefficients) <- coef.names
names(turnout.npopp$coefficients) <- coef.names
names(turnout.udgov$coefficients) <- coef.names
names(turnout.udopp$coefficients) <- coef.names
# Create joint table
turnout.other.part <- mtable(turnout.npgov, turnout.npopp, turnout.udgov, turnout.udopp,
														 summary.stats=c("N"))


